function fm = fm_with_neighborhood (c1, c2, rad)
	C = contingency_matrix(c1, c2, rad);	  %form contingency matrix	
	n = length(c1.bmus); %data quantity

	sumC=sum(sum(C.^2));	 %sum over rows & columnns of nij^2
	ni = sum(C,2);
	ni = ni.*(ni-1)/2;
	nis = sum(ni);
	nj = sum(C,1);
	nj = nj.*(nj-1)/2;
	njs = sum(nj);
	fm = 0.5*(sumC-n)/sqrt(nis*njs);    % FM - Fowlkes and Mallows 1983
	